home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 41 / Amiga Format CD41 (1999-06)(Future Publishing)(GB)[!][issue 1999-07].iso / -seriously_amiga- / programming / other / scm / slib / root.scm < prev    next >
Text File  |  1999-04-19  |  7KB  |  205 lines
  1. ;;;"root.scm" Newton's and Laguerre's methods for finding roots.
  2. ;Copyright (C) 1996, 1997 Aubrey Jaffer
  3. ;
  4. ;Permission to copy this software, to redistribute it, and to use it
  5. ;for any purpose is granted, subject to the following restrictions and
  6. ;understandings.
  7. ;
  8. ;1.  Any copy made of this software must include this copyright notice
  9. ;in full.
  10. ;
  11. ;2.  I have made no warrantee or representation that the operation of
  12. ;this software will be error-free, and I am under no obligation to
  13. ;provide any services, by way of maintenance, update, or otherwise.
  14. ;
  15. ;3.  In conjunction with products arising from the use of this
  16. ;material, there shall be no use of my name in any advertising,
  17. ;promotional, or sales literature without prior written consent in
  18. ;each case.
  19.  
  20. (require 'logical)
  21.  
  22. ;;;; Newton's Method explained in:
  23. ;;; D. E. Knuth, "The Art of Computer Programming", Vol 2 /
  24. ;;; Seminumerical Algorithms, Reading Massachusetts, Addison-Wesley
  25. ;;; Publishing Company, 2nd Edition, p. 510
  26.  
  27. (define (newton:find-integer-root f df/dx x_0)
  28.   (let loop ((x x_0) (fx (f x_0)))
  29.     (cond
  30.      ((zero? fx) x)
  31.      (else
  32.       (let ((df (df/dx x)))
  33.     (cond
  34.      ((zero? df) #f)        ; stuck at local min/max
  35.      (else
  36.       (let* ((delta (quotient (+ fx (quotient df 2)) df))
  37.          (next-x (cond ((not (zero? delta)) (- x delta))
  38.                    ((positive? fx) (- x 1))
  39.                    (else (- x -1))))
  40.          (next-fx (f next-x)))
  41.         (cond ((>= (abs next-fx) (abs fx)) x)
  42.           (else (loop next-x next-fx)))))))))))
  43.  
  44. (define (integer-sqrt y)
  45.   (newton:find-integer-root (lambda (x) (- (* x x) y))
  46.                 (lambda (x) (* 2 x))
  47.                 (ash 1 (quotient (integer-length y) 2))))
  48.  
  49. (define (newton:find-root f df/dx x_0 prec)
  50.   (if (and (negative? prec) (integer? prec))
  51.       (let loop ((x x_0) (fx (f x_0)) (count prec))
  52.     (cond ((zero? count) x)
  53.           (else (let ((df (df/dx x)))
  54.               (cond ((zero? df) #f) ; stuck at local min/max
  55.                 (else (let* ((next-x (- x (/ fx df)))
  56.                      (next-fx (f next-x)))
  57.                     (cond ((= next-x x) x)
  58.                       ((> (abs next-fx) (abs fx)) #f)
  59.                       (else (loop next-x next-fx
  60.                               (+ 1 count)))))))))))
  61.       (let loop ((x x_0) (fx (f x_0)))
  62.     (cond ((< (abs fx) prec) x)
  63.           (else (let ((df (df/dx x)))
  64.               (cond ((zero? df) #f) ; stuck at local min/max
  65.                 (else (let* ((next-x (- x (/ fx df)))
  66.                      (next-fx (f next-x)))
  67.                     (cond ((= next-x x) x)
  68.                       ((> (abs next-fx) (abs fx)) #f)
  69.                       (else (loop next-x next-fx))))))))))))
  70.  
  71. ;;; H. J. Orchard, "The Laguerre Method for Finding the Zeros of
  72. ;;; Polynomials", IEEE Transactions on Circuits and Systems, Vol. 36,
  73. ;;; No. 11, November 1989, pp 1377-1381.
  74.  
  75. (define (laguerre:find-root f df/dz ddf/dz^2 z_0 prec)
  76.   (if (and (negative? prec) (integer? prec))
  77.       (let loop ((z z_0) (fz (f z_0)) (count prec))
  78.     (cond ((zero? count) z)
  79.           (else
  80.            (let* ((df (df/dz z))
  81.               (ddf (ddf/dz^2 z))
  82.               (disc (sqrt (- (* df df) (* fz ddf)))))
  83.          (if (zero? disc)
  84.              #f
  85.              (let* ((next-z
  86.                  (- z (/ fz (if (negative? (+ (* (real-part df)
  87.                                  (real-part disc))
  88.                               (* (imag-part df)
  89.                                  (imag-part disc))))
  90.                         (- disc) disc))))
  91.                 (next-fz (f next-z)))
  92.                (cond ((>= (magnitude next-fz) (magnitude fz)) z)
  93.                  (else (loop next-z next-fz (+ 1 count))))))))))
  94.       (let loop ((z z_0) (fz (f z_0)) (delta-z #f))
  95.     (cond ((< (magnitude fz) prec) z)
  96.           (else
  97.            (let* ((df (df/dz z))
  98.               (ddf (ddf/dz^2 z))
  99.               (disc (sqrt (- (* df df) (* fz ddf)))))
  100.          ;;(print 'disc disc)
  101.          (if (zero? disc)
  102.              #f
  103.              (let* ((next-z
  104.                  (- z (/ fz (if (negative? (+ (* (real-part df)
  105.                                  (real-part disc))
  106.                               (* (imag-part df)
  107.                                  (imag-part disc))))
  108.                         (- disc) disc))))
  109.                 (next-delta-z (magnitude (- next-z z))))
  110.                ;;(print 'next-z next-z )
  111.                ;;(print '(f next-z) (f next-z))
  112.                ;;(print 'delta-z delta-z 'next-delta-z next-delta-z)
  113.                (cond ((zero? next-delta-z) z)
  114.                  ((and delta-z (>= next-delta-z delta-z)) z)
  115.                  (else
  116.                   (loop next-z (f next-z) next-delta-z)))))))))))
  117.  
  118. (define (laguerre:find-polynomial-root deg f df/dz ddf/dz^2 z_0 prec)
  119.   (if (and (negative? prec) (integer? prec))
  120.       (let loop ((z z_0) (fz (f z_0)) (count prec))
  121.     (cond ((zero? count) z)
  122.           (else
  123.            (let* ((df (df/dz z))
  124.               (ddf (ddf/dz^2 z))
  125.               (tmp (* (+ deg -1) df))
  126.               (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf))))
  127.               (df+sqrt-H (+ df sqrt-H))
  128.               (df-sqrt-H (- df sqrt-H))
  129.               (next-z
  130.                (- z (/ (* deg fz)
  131.                    (if (>= (magnitude df+sqrt-H)
  132.                        (magnitude df-sqrt-H))
  133.                    df+sqrt-H
  134.                    df-sqrt-H)))))
  135.          (loop next-z (f next-z) (+ 1 count))))))
  136.       (let loop ((z z_0) (fz (f z_0)))
  137.     (cond ((< (magnitude fz) prec) z)
  138.           (else
  139.            (let* ((df (df/dz z))
  140.               (ddf (ddf/dz^2 z))
  141.               (tmp (* (+ deg -1) df))
  142.               (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf))))
  143.               (df+sqrt-H (+ df sqrt-H))
  144.               (df-sqrt-H (- df sqrt-H))
  145.               (next-z
  146.                (- z (/ (* deg fz)
  147.                    (if (>= (magnitude df+sqrt-H)
  148.                        (magnitude df-sqrt-H))
  149.                    df+sqrt-H
  150.                    df-sqrt-H)))))
  151.          (loop next-z (f next-z))))))))
  152.  
  153. (define (secant:find-root-1 f x0 x1 prec must-bracket?)
  154.   (letrec
  155.       ((stop?
  156.     (cond ((procedure? prec) prec)
  157.           ((and (integer? prec) (negative? prec))
  158.            (lambda (x0 f0 x1 f1 count)
  159.          (>= count (- prec))))
  160.           (else
  161.            (lambda (x0 f0 x1 f1 count)
  162.          (< (min (abs f0) (abs f1)) prec)))))
  163.        (bracket-iter
  164.     (lambda (xlo flo xhi fhi count)
  165.       (if (stop? xlo flo xhi fhi count)
  166.           (if (> (abs flo) (abs fhi)) xhi xlo)
  167.           (let tail ((del (- (/ flo (- fhi flo)))))
  168.         (let* ((xnew (+ xlo (* del (- xhi xlo))))
  169.                (fnew (f xnew)))
  170.           (cond ((< 0.1 del 0.9)    ;accept secant step
  171.              (if (positive? fnew)
  172.                  (bracket-iter xlo flo xnew fnew (+ count 1))
  173.                  (bracket-iter xnew fnew xhi fhi (+ count 1))))
  174.             (else            ;bisection step
  175.              (tail 1/2)))))))))
  176.  
  177.     (let ((f0 (f x0))
  178.       (f1 (f x1)))
  179.       (cond ((<= f0 0 f1)
  180.          (bracket-iter x0 f0 x1 f1 0))
  181.         ((<= f1 0 f0)
  182.          (bracket-iter x1 f1 x0 f0 0))
  183.         (must-bracket? #f)
  184.         (else
  185.          (let secant-iter ((x0 x0)
  186.                    (f0 f0)
  187.                    (x1 x1)
  188.                    (f1 f1)
  189.                    (count 0))
  190.            (cond ((stop? x0 f0 x1 f1 count)
  191.               (if (> (abs f0) (abs f1)) x1 x0))
  192.              ((<= f0 0 f1)
  193.               (bracket-iter x0 f0 x1 f1 count))
  194.              ((>= f0 0 f1)
  195.               (bracket-iter x1 f1 x0 f0 count))
  196.              ((= f0 f1) #f)
  197.              (else
  198.               (let ((xnew (+ x0 (* (- (/ f0 (- f1 f0))) (- x1 x0)))))
  199.             (secant-iter xnew (f xnew) x0 f0 (+ count 1)))))))))))
  200.  
  201. (define (secant:find-root f x0 x1 prec)
  202.   (secant:find-root-1 f x0 x1 prec #f))
  203. (define (secant:find-bracketed-root f x0 x1 prec)
  204.   (secant:find-root-1 f x0 x1 prec #t))
  205.